# ------------------------------------- #
# Replication code for:
#
# Rathbun, Brian, Christopher Sebastian Parker, and Caleb Pomeroy "Separate but Unequal: Ethnocentrism and Racialization 
# Explain the 'Democratic' Peace in Public Opinion," American Political Science Review.
#
# This script reproduces our reanalyses of Michael Tomz and Jessica Weeks (2013) "Public Opinion and 
# the Democratic Peace," American Political Science Review.
# ------------------------------------- #

# --- libraries --- #
library(psych)
library(ggeffects)
library(ggplot2)
library(texreg)
library(ltm)
set.seed(1912)

# --- set working directory
setwd("~/Downloads/replication_files/")

# --- import data
tomz_weeks <- readRDS("data/tomz_weeks.rds")
nrow(tomz_weeks) # N = 1273, as mentioned in the paper's appendix 

# --- ethnocentrism measure
# reduce ethnocentrism items to unidimensional factor scores
ethno_data <- tomz_weeks[,c("immig_ls", "immig_bp", "immig_q", "affaxn")]
cronbach.alpha(na.omit(ethno_data)) # as footnoted, cronbach alpha = 0.678
ethno_fit <- fa(ethno_data, nfactors=1, rotate="varimax", scores=TRUE, fm="pa")
ethno_fit$loadings # as footnoted, SS loadings = 1.672
tomz_weeks$ethno_scale <- ethno_fit$scores[,1] 

# also include a simple additive version of ethnocentrism
tomz_weeks$ethno_scale_additive <- rowSums(tomz_weeks[,c("immig_ls", "immig_bp", "immig_q", "affaxn")]) 

# --- regression models
summary(mod_base <- glm(strikef ~ democ + ally + trade + male + white + age + college + hawk + ethno_scale,
                        data = tomz_weeks,
                        family = "binomial"))

summary(mod_ethno <- glm(strikef ~ democ + ally + trade + male + white + age + college + hawk + ethno_scale*democ,
                         data = tomz_weeks,
                         family = "binomial"))

summary(mod_hawk <- glm(strikef ~ democ + ally + trade + male + white + age + college + ethno_scale + hawk*democ,
                        data = tomz_weeks,
                        family = "binomial"))

# --- table A2
screenreg(list(mod_base, mod_ethno, mod_hawk))

# as footnoted in the appendix, substantively similar interaction between additive ethnocentrism and regime type treatment,
# albeit the p-value increases to 0.057
summary(glm(strikef ~ democ + ally + trade + male + white + age + college + hawk + ethno_scale_additive*democ,
            data = tomz_weeks,
            family = "binomial"))

# --- ethnocentrism x regime type interaction plot
ethno_predict <- ggpredict(mod_ethno, c("ethno_scale", "democ")) 

plot_df <- data.frame(strike = ethno_predict$predicted,
                      ethno = ethno_predict$x,
                      democ = ethno_predict$group,
                      ci_high = ethno_predict$conf.high,
                      ci_low = ethno_predict$conf.low)

# as mentioned in the appendix, 18.1% max vs min difference in support for strike 
max(subset(plot_df, democ == 1)$strike)-max(subset(plot_df, democ == 0)$strike)

# --- figure 4 and A4
ggplot(plot_df, aes(x = ethno, y = strike, group = democ)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .1) +
  geom_line(aes(color = democ), size = 2) +
  scale_color_manual(labels = c("No", "Yes"), values=c("#92978a","gray30",  "#c5bec2")) +
  labs(y = "Support for Strike", x = "Ethnocentrism", color = "Democracy?") +
  theme_minimal() +
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 18, margin = margin(t=10)),
        axis.title.y = element_text(size = 18, margin = margin(r=10)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        legend.key = element_rect(color = NA, fill = NA),
        legend.key.size = unit(1, "cm"))

# --- conditional average treatment effects
tomz_weeks$ethno_scale_binary <- ifelse(tomz_weeks$ethno_scale > median(na.omit(tomz_weeks$ethno_scale)), 1, 0)
# the democ coef represents respondents at or below the ethno median;
# the democ:ethno_scale_binary coef represents respondents above the ethno median
summary(mod_ethno_cate <- glm(strikef ~ democ + ally + trade + male + white + age + college + hawk + ethno_scale_binary*democ,
                         data = tomz_weeks,
                         family = "binomial"))

# --- panel analysis, i.e. within subjects results
tomz_weeks_panel <- subset(tomz_weeks, post == 1)
# assign a 1 to subjects who did not support strikes against democracies but supported strikes against nondemocracies, 0 otherwise
tomz_weeks_panel$dempeace_panel <- ifelse(tomz_weeks_panel$democ == 1 & tomz_weeks_panel$strikef == 0 & tomz_weeks_panel$democ_post == 0 & tomz_weeks_panel$strikef_post == 1 |
                                            tomz_weeks_panel$democ == 0 & tomz_weeks_panel$strikef == 1 & tomz_weeks_panel$democ_post == 1 & tomz_weeks_panel$strikef_post == 0,
                                    1, 0)

summary(mod_ethno_panel <- glm(dempeace_panel ~ ally + trade + male + age + ethno_scale, 
                               family = "binomial", 
                               data = tomz_weeks_panel))
# --- table A3
screenreg(mod_ethno_panel)

# estimate predicted values
ethno_predict_panel <- ggpredict(mod_ethno_panel, "ethno_scale")

plot_df <- data.frame(strike = ethno_predict_panel$predicted,
                      ethno = ethno_predict_panel$x,
                      democ = ethno_predict_panel$group,
                      ci_high = ethno_predict_panel$conf.high,
                      ci_low = ethno_predict_panel$conf.low)
max(plot_df$strike) # nearly 2.5 times larger, as mentioned in the appendix text
min(plot_df$strike)

# --- figure A5
ggplot(plot_df, aes(x = ethno, y = strike)) +
  geom_ribbon(aes(ymax = ci_high, ymin = ci_low), alpha = .2, fill = "black") +
  geom_line(color ="gray30" , size = 2) +
  labs(y = "Predicted Democratic\nPeace Tendencies", x = "Ethnocentrism") +
  theme_minimal() +
  theme(legend.title = element_text(size = 16),
        legend.text = element_text(size = 16),
        axis.text.y = element_text(size = 16),
        axis.text.x = element_text(size = 16),
        axis.line.y.right = element_blank(),
        axis.title.x = element_text(size = 18, margin = margin(t=10)),
        axis.title.y = element_text(size = 18, margin = margin(r=10)),
        panel.spacing.x = unit(1.5, "lines"),
        panel.grid.minor.x = element_blank(),
        legend.key = element_rect(color = NA, fill = NA),
        legend.key.size = unit(1, "cm"))

